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We theoretically study the dynamic polarization of lattice nuclear spins in GaAs double quantum 
dots containing two electrons. We introduce a semiclassical model that allows to explore a wide range 
of parameter regimes in this system. We identify three regimes of long-term dynamics, including the 
build up of a large difference in the Overhauser fields across the dots, the saturation of the nuclear 
polarization process associated with formation of so-called "dark states," and the elimination of the 
difference field. When the dots are different sizes, build up of difference fields follows the polarization 
process, whereas for nearly identical dots, we find the dynamics are strongly influenced by nuclear 
spin noise. In the absence of noise, all three steady states are achieved depending on parameters. 
However, with nuclear spin noise we find the build up of difference fields competes with polarization 
saturation in dark states, while the elimination of the difference field does not, in general, correspond 
to a stable steady state of the polarization process. These results are in agreement with dynamic 
nuclear polarization experiments in double quantum dots. 
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I. INTRODUCTION 

The study of non-equilibrium dynamics of nuclei in 
solids has a long histor}^ and has become particularly 
relevant as nanoscale engineering and improvements in 
control allow to probe mesoscopic collections of nuclear 
spins 2 8 . This control has direct applicability to quantum 
information science, where nuclear spins are often a main 
source of dephasing^. The goal of developing an under- 
standing of electronic control of nuclei is to circumvent 
this nuclear dephasing and to turn nuclear spins into a 
useful resource 10 , as indicated in recent experiments^^. 

Double quantum dots in III-V semiconductors can be 
operated with two electrons coupled to approximately 
10 4 to 10 6 nuclei by the contact hyperfine interaction. 
Repeated cycles transitioning from the electronic singlet 
to triplet states can be used to polarize the nuclear spins; 
electron spin flips between the singlet and triplet spaces 
occur due to the difference D in the Overhauser fields on 
the two dots^. Early experimental and theoreticaJ^HUl 
work suggested that the polarization process naturally 
drove the projection of the difference field onto the mag- 
netic field axis D z to zero. However, later experiments 
and theory both showed that the polarization is naturally 
accompanied by a growth in D z and that the data in the 
original experiment s show ing a suppression in D z was 
likely misinterpreted^ 2 ^. Instead the results are more 
consistent with the growth of a large D z accompanied by 
a reduction in measurement contrast between singlet and 
triplet states, which makes it appear as if D z is smalP^. 

In Ref. 23 we developed a model to describe the long 
time dynamics of the nuclear spins undergoing adiabatic 
pumping. These results are in good agreement with the 
experiments described above. In the present work we 
present a detailed theoretical analysis of these problems. 
We describe the theoretical methods developed to study 



this system and give a systematic description of the po- 
larization phenomena relevant to experiments. 

Our theoretical methods are based on a semiclassical 
description of the nuclear spin dynamics in which the 
nuclear spins are grouped into small sets, each homoge- 
neously coupled to the electron spirP^. The nuclei in each 
set may be treated as a single collective spin and a semi- 
classical treatment is justified provided the number of 
spins in each set remains large. Increasing the number of 
such sets improves the approximation to the true hyper- 
fine coupling. More formally, we construct a systematic 
approximation to the true hyperfine coupling in terms of 
a reduced set of M coupling constants. For the optimal 
choice of coupling constants, we rigorously prove that our 
approximation reproduces the exact semiclassical time 
dynamics to within a fixed error for a time that increases 
linearly with M . For large M . this allows examination of 
the long timescales relevant for polarization experiments. 
This approach extends previous work that assumes that 
all nuclei o n a gi ven dot have equal coupling to the elec- 
tron S pi n lM22l26H2E an a pp r oach which often incorrectly 
predicts rapid saturation of the polarization. Other ex- 
tensions to this homogenous coupling m odel, including 
semiclassical solutions for the central spir l 29 * 30 *, and clus- 
ter and diagramatic expansion techniques for short time 
non-equilibrium behavioi^^ do not explore the wide 
range of time scales or relevant physics for the double 
dot case. 

Our results can be broken up into two distinct cases de- 
pending on whether or not the dots are identical. When 
the dots are different sizes, then the hyperfine coupling, 
which scales inversely with the volume, is larger on the 
smaller dot and we find that the Overhauser field grows 
preferentially on the smaller dot as the polarization in- 
creases. This preferential growth results in a large Over- 
hauser difference field D z . For two dots with a difference 
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in volume of less than ~ 20% we find a rich and com- 
plex phase diagram for the nuclear spin dynamics, which 
can be broken into two distinct regimes. The first regime 
occurs with large external magnetic fields or short cycle 
times. In this regime the system saturates without signif- 
icant polarization because the perpendicular components 
of D rapidly approach zero and spin flips are suppressed; 
the system approaches a semiclassical "dark state." This 
occurs with no statistical change in the distribution of 
D z . The second regime occurs in the limit of smaller 
magnetic fields or slower cycle times. In this regime, the 
dynamics are sensitive to the inclusion of nuclear spin 
noise. In the absence of nuclear spin noise we find one 
potential end state of polarization is a "zero state" in 
which all components of D — >• 0. In this state the sin- 
glet and triplet electronic subspaces are completely de- 
coupled and spin flips no longer occur. Simultaneously, 
though, there are instabilities leading to the growth of 
large Overhauser difference fields. Crucially, when even 
a small amount of nuclear spin noise is added the zero 
states strongly destabilize and the system generically be- 
comes unstable to the growth of large difference fields as 
shown in Ref. [23] 

These results provide a clear picture of the polarization 
dynamics in such double quantum dot systems and will 
be a useful guide to future experiments aimed at more 
precise control of the nuclear spins. Although the paper 
is specific to double quantum dots in GaAs, many of the 
results and theoretical methods extend to other central 
spin systems under invest igatiorP^SD. More generally, 
this work is of fundamental interest as we explore the 
dynamics of an interacting, many-body system when it 
is far from equilibrium 37 . 

The paper is organized as follows. In section II we 
define the Hamiltonian for the double dot system and in- 
troduce the polarization cycle. In section III we system- 
atically derive a semiclassical model for the nuclear spins 
starting from the coarse-grained evolution of the nuclear 
spin density matrix. In section IV we present our results 
for identical and unequal dots in the presence and ab- 
sence of nuclear spin noise. In appendix A we provide a 
summary of the parameters used in our simulations. In 
appendix B we describe our approach to coarse graining 
the electron wave function and provide rigorous bounds 
on the error in time evolution due to the course graining. 
In appendix C we extend our simulations to the case of 
multiple nuclear species and find qualitatively the same 
results as for a single species. 



II. SETUP 

For a double quantum dot with two electrons, we can 
write the Hamiltonian for the lowest energy (1,1) and 
(0, 2) electron states, where (n, m) indicates n (m) elec- 
trons in the left (right) dot. To model nuclear polariza- 
tion, we first derive an effective two-level Hamiltonian 
to describe the system near the crossing of the singlet s 




FIG. 1: a) The Overhauser field in each dot gives rise to sum 
and difference fields which are relevant for the double dot sys- 
tem, b) Schematic of two-electron energy levels as a function 
of detuning e between (1,1) and (0,2) charge states. Arrows 
indicate adiabatic sweep through avoided crossing (pink) and 
rapid sweep back to (0,2) with reload (green), c) Spin-flip 
pathways between the s and T + states as the exchange en- 
ergy J(e) is swept through the crossing, showing the nuclear 
operators involved in each path. Each pathway is a term in 
D- in Eq.[5] 

and lowest energy triplet state, T + , of this two-electron 
system, then solve the time dynamics. Dynamic nuclear 
polarization (DNP) experiments operate near this cross- 
ing, typically with an adiabatic sweep of the difference in 
the dots electric potential through the s-T + degeneracy 
(Fig.[]Ji), followed by a non-adiabatic return to (0,2) and 
reset of the electronic state via coupling to leads. 

If ^(r) is the single-particle envelope wave function 
on dot d = Z,r (for the left, right dot), the effective 
hyperfine coupling for the nuclear spin at r^d is 9kd = 
^hf^\^d{^kd)\ 2 where a,hf is the hyperfine coupling con- 
stant, and vo is the volume per nuclear spin. We in- 
troduce two collective nuclear spin operators to denote 
the Overhauser fields in the left (L) and right (R) dots, 
L = llkdki^ki and R = Y^kdkAkr, and further define 
S = (L + R)/2, D = (L - R)/2, where I kd is the an- 
gular momentum of the k th nucleus on dot d. The rms 
Overhauser energy in the infinite temperature ensemble 
is Qd — (52k9kdI(I + l)/^) 1 ^ 2 where / is the magnitude 
of each nuclear spin. We define ft = (0| + Q!*)/2, and 
choose units such that Q = — g ^ B = 1, where g* is the 
electron effective g-factor and \±b is the Bohr magneton. 
In the basis {\s) , |T+> , |T ) , |T_», where the T m are the 
(1,1) triplet states and s is the (l,l)-(0,2) hybridized 
singlet state, the Hamiltonian is^ 

/ -J(e) vt) + -V2vD z -vD- \ 

H= vt)_ -B ext + S z S-/V2 

-V2vD z S+/V2 S-/V2 

\ -vt) + £ext - Sj 

where D± = D x ±iD y and similarly for S±, B ext is an ex- 
ternal magnetic field, and v = v (e) = cos 0(e)/ y/2. The 
parameters cos 9(e), the overlap of the adiabatic singlet 
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state with the (1,1) singlet state, and J(e), the splitting 
between s and To, are both functions of the energy differ- 
ence s between the (1,1) and (0,2) charge states. Here 
the nuclear spin variables refer to the full quantum me- 
chanical operators on the nuclear spin space. In appendix 
C we will consider the case of multiple nuclear species, 
but for now we consider the nuclei to be spin-3/2 of a 
single species, in a frame rotating at the nuclear Larmor 
frequency. 

Assuming that J, B ext > (], we perform a formal ex- 
pansion in the inverse electron Zeeman energy operator 
m = Q/(B ext — S z + if]). We apply a unitary transfor- 
mation that rotates the quantization axis of the triplet 
states to align with B ext — S and find the Hamiltonian 
for the {\s} , |T+)} subspace to second order in J, m: 



eff 



-J(e) + h a v(e)D + \ m 
v(e)D_ -B ext + hr)' K ) 



where 



9?; 2 ?; 2 

hs = —rD\D z ~ D- —— ¥ D + , (2) 

fiT — S z — -(S-S+m + mS-S+), 

b_ =D_ +mS-D z - \mS 2 mf) + - \mS-S+mD-, 
4 ^4 ^ 

D Z =D Z - ^(S+mD_ + S-mf) + ). 

Of particular interest is that the off-diagonal term, which 
produces nuclear polarization, vanishes in the semiclas- 
sical limit of (D) — >> 0, i.e., in the zero states. 



III. MODEL 

We develop a model for the evolution of the nuclear 
spin density matrix after one pair of electrons has passed 
through the system. We approximate the sweep through 
the \s)-\T+) degeneracy as a Landau-Zener process, which 
we solve approximately for the effect on the nuclear sys- 
tem. By coarse-graining this evolution over a cycle we de- 
rive a master equation for the nuclear spins. Finally, we 
add the effects of nuclear dipole-dipole interactions and 
quadrupole splittings phenomenologically. The deriva- 
tion presented here is complementary to that of Ref. [23] 
and results in the same equations of motion. 

The electron system is prepared in \s) at large negative 
t = — T/2, where T is the total cycle time. We identify 
the (nuclear spin) eigenstates of the operator D+D_, la- 
beled \D±) with eigenvalues D\. For initial state |^o) = 
\s) (g) \D±), the crossing either leaves the electronic state 
unchanged or flips an electron and nuclear spin to the 
state |T+) <g> \D f ± ), where \D' ± ) = D^D- \D±). We note 
that \D'j_) is an eigenstate of D-D+ with eigenvalue D\. 
The problem is now reduced to finding Landau-Zener so- 
lutions for each independent two-level system {\s)(g>\D±), 



|T + ) <S> \D' ± )}, where we approximate h s — » (D± \ h s \D±), 
Jit —> (D' ± \ fir \D' ± ) : valid to leading order in 1/J and 
1/N. We model the actual sweep of e by a linear sweep of 

J so J(t) = -2/3 2 £ + £ ext , where /? = ^\dJ(e)/dt\\ t = . 

We take v(e) « 1/V2 to be constant, valid in the limit 
of large tunnel coupling, and assume j3 <C B ext to ensure 
the applicability of Eq. [I] In this limit, the off-diagonal 
part of iY e ff in Eq. [I] produces standard Landau-Zener 
behavior, while the diagonal components of H e R are sim- 
ply phases picked up by the nuclei, depending on which 
electronic state is occupied. 

After one cycle, |^o) evolves into \^\) = cs \s)®\D±_) + 
ct \T+)®\D'±). For f3 2 T > 1, the standard Landau-Zener 
formula gives the flip probability asp/ = 1— exp(— 27ru; 2 ), 
where uo = v(D±)//3, and 



cs = V 1 ~Pf exp(-z0 s ), c T 
(j) S « / h s dt 



^7 exp(-z0 T ) 



T/2 
o 



(3) 



/to 
-T/2 



h s dt + (T/2 - t )h T + 0adM, 



where the crossing occurs at a time to ~ S z / ' f3 2 . We in- 
clude in <pT the phase picked up by following the adiabat, 
4>ad- We approximate <\>ad by interpolating between the 
limits uj = v(D±)//3 —> and ou — >• oo, giving^ 



: 2ttu/ 



Pf 



1 - 2tt + log 



where r = T/3/2. 

We move from the independent two-level systems to 
the general case by noting that the components of \Sfr) 
depend only on the eigenvalue D± and on the polariza- 
tion S z (which we approximate as commuting). Thus 
we can write an operator pf = ^2r> ± Pf(D_\_) \D±) (D±\, 

and similarly for 4>Si4>t- The nuclear spin density ma- 
trix after each cycle is given by tracing over the electronic 
states. The nuclear density matrix evolution is then 



V 1 ~Pf 



Pf 



D+D. 



-D 



+ h 



where p n is the nuclear density matrix after n cycles. 

Rather than solve for the exact dynamics of the nuclear 
density matrix-still an intractably hard computational 
problem for any reasonable number of nuclear spins-we 
instead adopt an approximate solution to the problem 
using the P-representation for the density matrix as an 
integral over products of spin coherent states. From the 
thermal distribution, we choose such a spin coherent state 
and evolve it, where we interpret expectation values (...) 
as being taken in that state. The ensemble of such tra- 
jectories represents the physical system^. 



4 



We organize this calculation by noting that the com- 
ponents of the Landau-Zener model ((t>Si<t>T,Pf,D±) are 
only functions of L and R. A spin coherent state is en- 
tirely described by its expectation values = (1^). For 
the kth spin on the left dot, we expand the discrete time 
difference (Iki) n — (Ifcz) n _i after n and n—1 cycles in the 
small parameter g k i, giving an evolution equation 



dhi 
dt 



= Qki^Pi^ (i[d gkl L^l k i)^j = g M Pi x \ kh (4) 

/x=l 



with 



1 

f 



(l-Pf)(Vi<l>s) + (p»(V/0 T > - Im(7i) 
where Vi = (d Lx ,d Ly ,d Lz ) and 



Pf 



D_D + 



(5) 



and similarly for i/~ r , P r , and 7 r , with L replaced by R. 
The factorization of expectation values is a natural con- 
sequence of our spin-coherent state approximation, as it 
explicitly prevents entanglement between spins. Thus we 
have an effective, semi-classical picture of nuclear spins 
precessing and being polarized by their interaction with 
the electron spin, integrated over one cycle. 

We approximate the electron wavefunction as a 
piecewise-flat function with M levels, which we refer to 
as the annular approximation, as illustrated in Fig. [2^l. 
Each annulus defines l n d = ^Z ken hd^ where the sum is 
over all nuclei with the same hyperfine coupling to the 
electron. Since g k is identical for all k G n, we can sim- 
ply replace i k d with I n d in Eq. [4] Furthermore, I 2 is 
a conserved quantity, so we can study the evolution of 
M <C N spins in a red uced Hilbert space. The typical 
size of I n is ~ y^N/M ^> 1, which allows us to replace 
the spin-coherent states used above with semi-classical 
spins, and makes taking expectation values straightfor- 
ward: all quantum operators can be replaced by their 
expectation values directly. The annular approximation 
should correctly describe the nuclear dynamics for a time 
scale given by the inverse of the difference between the 
g k of adjacent annuli. 

To illustrate, to first order in mo = B^., for d = Z, r, 



P d =p f X (A+i - A S ± ) + moT y^z x D (6) 

f3 2 

^RPf^d^ADT 1 



Pf D z 
27TLJ 2 

Im -7 r ) 



4ttv 2 

+ (1 - p f X/2)(A D z z + A_D ± ) 

where the top sign applies for d — = (D x , D y ,0), 

S ± = (S x ,S y ,0), A = 1 - 2t /T gives the shift in the 
location of the crossing, and Ao, A_, A + , Ao, and 
Tq are constants depending on the details of the pulse 
cycle (see below). We have replaced operators by their 




Ao — |r -» 
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L±l \T + ) 
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A+ — |T -» 
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r - - lT - } 
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\ T +) 
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A — ,T - } 
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\T + ) 

\S) 
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FIG. 2: a) Independent Random Variable Annular Approx- 
imation (IRVAA) to the electron wavefunction in the double 
dot. b) Key processes contributing to Eq. [6] 



expectation values and removed the angle brackets since 
we are now in the semiclassical limit. To leading order 
in mo, Im(7i — 7 r ) = 2(D x z)pf /D\. It is clear from 
Eqs. [4][6] that all dynamics stop in the zero states with 
D = 0, consistent with the idea that true saturation of 
polarization requires that all components of D be small. 
We will focus on the stability of such states in various 
parameter regimes. The equations of motion in Ref. [23] 
are found from Eq. [6] by including only the lowest order 
in Q T and Q//3, which is the limit of fast cycles and small 
spin flip probability per cycle, respectively. 

First we outline the meanings of the parameters in the 
model. As indicated schematically in Fig.[2]3, the Tq term 
originates in the hyperfine flip-flop, the Ao and A_ terms 
are the off-resonant effects of coupling from the singlet 
state to the To and T_ states, respectively, Ao comes from 
coupling between the T + and Tq states, and A + comes 
from Knight shifts due to occupation of the T + state. To 
leading order in mo, for a pulse sequence consisting of 
only the Landau-Zener sweep, with instantaneous eject 
and reload, the parameters have values 



m /4 



A = 




a ra , A_ 




W)/c 


\j(t) 


A+ = 


1/4, 


Ao 


= m /4 


r = 


p 2 ' 




= fc 



where f c = 1/T is the cycle frequency and (.) c indicates 
an average taken over a full cycle; these values can be 
modified readily by changing the details of the pulse cy- 
cle, while leaving the Landau-Zener portion unchanged. 
In Appendix A we provide a reference for all parameters 
used in the simulations. 

Equation [4] is a good approximation of the nuclear dy- 
namics over a few DNP cycles because other nuclear pro- 
cesses are slow compared to a typical experimental cycle 
(~ 10- 100 ns^J). However, the full DNP may last millions 
of cycles at which point these other nuclear processes be- 
come important. Apart from Larmor precession, which 
is only relevant for the case of multiple nuclear species 
considered in Appendix C, nuclear quadrupole splittings 
and nuclear dipole-dipole interactions are the dominant 
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processes. They become relevant on a timescale of a few 
hundred microseconds in these system^. We include 
them in our model phenomenologically by adding a fluc- 
tuating magnetic field hkd(t) in the z-direction at each 
site (the transverse terms are strongly suppressed by the 
external field), such that 

= gtd-Pd x i k d ~ 7n hkd z x i k d (7) 

where j n is the nuclear gyromagnetic ratio. We further 
assume that the this field can be treated as noise and 
characterized by a Gaussian, uncorrelated white noise 
spectrum 

ll(Kd(t)K>d'(t')) n = 2 V S(t - t')5 kk ,5 dd , (8) 
where (-) n are averages over the noised 

IV. RESULTS 

The polarization dynamics display three characteristic 
behaviors: growth of large difference fields, saturation in 
nuclear dark states defined by D± = 0, and preparation 
in zero states D = which are global fixed points of the 
nuclear dynamics in the absence of noise. In Ref. [23] this 
system was studied in a restricted model focusing on the 
case where noise was present. Therein it was found that 
when the two dots have different hyperfine couplings the 
system generically grows large difference fields, while for 
identical dots, depending on parameters, the system is 
either unstable to the growth of large difference fields or 
saturates in dark states; however, the zero states were 
not found to be a relevant steady state in any parameter 
regime. In the present work we focus on extending the re- 
sults of Ref. [23] to a larger, more experimentally relevant, 
parameter regime by using equations of motion correct 
to second order in mo with a more complete model of the 
Landau-Zener sweep as described in the previous section. 
In addition, we consider the nuclear dynamics in the ab- 
sence of noise. We also present the full analytical calcu- 
lations which were omitted from Ref. [23j In all physical 
parameter regimes we find qualitatively consistent results 
with Ref. [23] however, for a limited, unphysical parame- 
ter regime we do find solutions to the equations of motion 
in the absence of noise where the zero state is uniformly 
reached starting from a completely uncorrelated nuclear 
spin ensemble. 

The simulations shown below were performed with the 
equations of motion correct to second order in mo with 
ipd( r ) a 2D Gaussian. Taking v 2 « 1/2, we estimate 
that for experiments performed with £? ext = 10 mT with 
T = 25 ns 11 , m « 0.18, T « 0.20, but the A and A 
terms depend on the rest of the cycle. In each of the 
simulations, we choose initial magnitudes and directions 
of the spins I n by a procedure equivalent to choosing 
initial directions for each of the N n spin-3/2 nuclei in 
the n th annulus and evaluating I n = ^2 ken ifc explicitly 
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FIG. 3: Phase diagram for the simplified model presented 
in Ref. 1231 At each value of parameters, twenty runs were 
started with D z = —2, S z = —10, and all other components 
chosen randomly according to the infinite temperature ensem- 
ble. The colorscale indicates how many of those runs ended 
with \D Z \ increased. The dark region is of saturation and 
the light region is of instability. The dashed line shows the 
prediction of the simple model of Eq. |30| which captures the 
phase boundary, especially at low A_/Ao- For parameters 
used, see Table I. 




lo Sio m o 

FIG. 4: Phase diagram as in Fig. [3] except with varying ex- 
ternal magnetic field and without any noise added. The pa- 
rameters were scaled with mo as shown in Section III. There 
is a clear boundary between saturation at large To and insta- 
bility at lower values of To , with appropriately large values of 
Ao and A_. See Table I for parameters. The symbols 'x' and 
'o' mark the parameters used for Fig. [7] below. 

(see Appendix B). The relationship between simulation 
time and laboratory time depends on the details of the 
pulse cycle, including pauses and reloads not considered 
explicitly here, but simulation time is roughly in units of 
g~a X , where g max « 2Q 2 /a hf is the largest value of g k , 
so t = 400 is approximately 10 ms. 

To organize our results we recall the phase diagram for 
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identical dots and the simplified model derived in Ref . [23] 
in which the only non-zero parameters are Ao,-, To and 
77, which corresponds to the limit of large magnetic fields 
and fast sweeps including nuclear spin noise. To obtain 
the phase diagram we consider for each set of parame- 
ters whether the system supports self-consistent growth 
of \D Z \ starting from large values of \D Z \ and \S Z \. This 
approach avoids complications with the metastability of 
zero states discussed later. Such simulations produce 
the phase diagram in Fig. 2b of Ref. 23, which is re- 
produced in Fig. 3 with the full data presented. From 
this figure it is clear that we can separate the dynamics 
into two regimes depending on parameters. For large ra- 
tios of To/Ao, which corresponds to large magnetic fields 
or strong pumping the system quickly saturates with no 
growth of large difference fields. For small ratios there is 
an instability towards large difference fields. In the first 
section we explore the dynamics in the absence of noise 
for identical dots with all parameters included. In the 
second section we include nuclear spin noise and asym- 
metry in the dot sizes. 



A. Noise Free Nuclear Spins 

From the general arguments given in the introduction 
it is clear that when the dots have different hyperfine cou- 
plings the system naturally grows a large difference field. 
Furthermore, in Ref. |23|it was shown that even identical 
dots display similar behavior in the presence of noise. In 
this section we analyze the case of identical dots in the 
absence of noise to better understand the role of the co- 
herent nuclear dynamics. We begin by deriving a phase 
diagram analogous to the one obtained in the presence 
of noise except we now look in the space of the experi- 
mentally accessible parameters cycle rate f c and inverse 
magnetic field mo- The results are shown in Fig. [4] where 
we see the same qualitative behavior as shown in Fig. [3] 
However, the dynamics are much richer than indicated by 
this simple phase diagram. In the following subsections 
we give examples of what happens to a nuclear spin en- 
semble starting from equilibrium for different parameters 
and regions of the phase diagram. 

Before proceeding, however, we note that in the ab- 
sence of noise the inhomogeneity of the electron wave- 
function plays a crucial role. This is because weak inho- 
mogeneity is equivalent to choosing the number of annuli 
M to be small and in this case the system moves rapidly 
to its maximally polarized state, with I n « — I n z for all 
n. Dynamics completely cease in this state, as can clearly 
be seen from Eq. [4j despite the fact that this state does 
not correspond to all of the nuclei being polarized, which 
would also require I n = 3N n /2. On the other hand, for 
strong inhomogeneity, or large M, when the system is 
not fully polarized other terms in compete with the 
polarization saturation and sustain the dynamics.^ 



1. Polarization Saturation 

When the magnetic field is large or the cycle rate is 
fast (i.e., Ao <C To), the system rapidly moves toward 
dark states (i.e., states with D± = 0), sending pf — >> 
without statistical change in the distribution of D Z1 as 
shown in Fig. [5^i. This limit is additionally characterized 
by only a small change in nuclear polarization as seen 
in Fig. [5)3. When the effects of the |s)-|Tq) coupling are 
important (i.e., Ao ~ To), the Ao term in Eq. [6] causes 
D± to increase, "rebrightening" the D± « dark states 
and allowing dynamics to continue. Coupling from the 
singlet to the Tq state is an essential ingredient in all 
of the effects discussed below. When Ao is significant, 
dynamics only stop near zero states with D = 0. 

2. Growth of difference fields 

Second, we observe the growth of large Overhauser 
fields. We consider a prototypical pulse sequence mo- 
tivated by experiments with moderate/large magnetic 
field, mo = 0.01 In this case, over 95% of the trajectories 
display a growth in \D Z \, as shown in Fig. [5^. We observe 
this behavior over a range of experimentally accessible 
magnetic fields and cycle frequencies. This increase in 
\D Z \ indicates that the spin flips are occurring predomi- 
nantly in one dot. We interpret these results as showing 
a continuing increase of \D Z \, where the peak of \D z (t) \ is 
an artifact of the annular approximation. Near the peak, 
many of the annular spins artificially reach their maximal 
polarization, at which point they should be broken into 
more annuli. Similar trajectories with different M show 
the maximum value of \D Z \ increasing with M (Fig.[5]i). 
The physical cause of this increase in \D Z \ is not clear, 
but it is associated with both A /r and A + /r being 
sufficiently large. When nuclear spin noise is included, 
the growth in (\D z \) e continues^. This could be the 
same phenomenon as seen in Ref. [T3l with transverse de- 
phasing helping to produce the large \D Z \ « B ext of that 
work, though unequal dot sizes could also produce that 
effectP. 



3. Zero states 

For moderate to small magnetic fields, when Ao ~ To, 
two different characteristic behaviors of particular note 
are observed. First, in the physical parameter regimes, 
which do not display general motion to zero states, the 
zero states are still important for the dynamics as they 
are a metastable state. That is, many trajectories spend 
a long time with \D Z \ near zero before escaping away to 
large \D Z \. This phenomenon is shown in the individual 
trajectory (thin red line) of Fig. [5^. 

Second, for parameters in our model which are not ex- 
perimentally accessible there is a mechanism that gives 
rise to attraction towards zero states. This is illustrated 
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FIG. 5: a-b) Simulations corresponding to the saturation region of the phase diagram. The solid lines are the median values of 
\D Z \ (a), S z (a) and D± (b) at each time step in an ensemble of 1000 trajectories. In all plots shaded regions show the 84 th and 
16 th percentiles, c) Simulations showing growth of (\D Z \) with the time shifted for each trajectory so that its maximum \D Z \ 
occurs at time zero. Bottom shows the median value of S z (S z ) e at each time step in an ensemble of 1000 trajectories. In the 
middle is similar (\D z \) e . Thin red line is a single trajectory. The curve at top shows the fraction of trajectories contributing 
to the ensemble at each time; this increases with time because some trajectories reach their maximum D z much later than 
others while the simulation time is fixed for each trajectory. 4.5% of the trajectories, which do not show this peak in \D Z \, are 
not included. Approximately 10% of the trajectories show behavior similar to that shown in the thin red line, where \D Z \ is 
reduced initially and then goes unstable to large \D Z \. d) Mean of the maximum value of \D Z \ reached on each trajectory for 
the same parameters as in (c) (open circles) except M varied between 20 and 160, with 5000 trajectories per point. Closed 
circles show similar results with mo = 0.05, r = 4 and all other parameters scaled appropriately. The physical system has 
M —> N « 10 6 , so we interpret this as an instability to large \D Z \, which is supported by simulations including transverse noise 
e) With different parameters, simulations showing reduction of (\D Z \) plotted as in (c) without the time 

D z 



IVB2) 



(see section 

shift. For these parameters, the trajectories have 



quickly, without time for strong polarization. 



in Fig. [5^, where we show an ensemble of trajectories in 
which D rapidly reduces toward zero. For the parameters 
of Fig. [5^, the standard deviation of D z was reduced by a 
factor of 28. We remark that as D — >• 0, the singlet state 
ceases mixing with the triplets and nuclear spin dynamics 
stop. Until something (outside this model, such as nu- 
clear dipole-dipole coupling) restores D, the polarization 
process is shut off, limiting the total nuclear polarization 
that can build up. While not shown in Fig. [5^, we observe 
a dramatic reduction of the total \D\, not just D Z1 consis- 
tent with this qualitative observation. However, because 
we have not observed this phenomenon in any physical 
parameter regimes we shall not study it further. 



4- Crossover 

For many choices of parameters, we find both trajecto- 
ries in which D z — >• and \D Z \ remains large, depending 
on initial conditions, as shown in Fig. |6^i. Note that 
when we add a small amount of transverse dephasing 
to these trajectories, as shown in Fig. [6}d, the median 
value of \D Z \ does not markedly change, but there are 



no longer trajectories with D z — >• 0; the noise apparently 
disrupts the fragile attraction toward \D Z \ —> 0. Simu- 
lations performed with parameters intended to approxi- 
mate experiment J 11 * 13 * are in this crossover regime. 

5. Stability of zero states 

We now investigate more carefully the stability of the 
zero states. Near the zero state the EOM are greatly 
simplified because many of the terms in arise from 
perturbative processes involving multiple applications of 
D. Keeping only the terms linear in D and working to 
first order in mo we can write 

D+ = (T + iA.)S*D + + (r moS*S+ - iA S* + )D z 

(9) 

D z = -Re[(r + iA_)D+S*_ ] - T m S ± • S* ± D Z , 

(10) 

where we have introduced the variable S* = 
^ kd gl d Ikd/2' Because dS / dt.dS* / dt ~ O(D), we can 
neglect the time dependence of S and in the EOM 
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1000 



1000 



FIG. 6: a) 1000 trajectories were run with initial conditions 
chosen from the thermal distribution with no noise. The mean 
value of \D Z \ is shown in black, and the the gray region en- 
closing 67% of the trajectories. A single trajectory is shown 
in the thin red line. For parameters, see Table 1. These pa- 
rameters are not represented in the phase diagram since they 
have very large A+. For these parameters, many trajectories 
are attracted near D = 0, as in the single trajectory shown, 
for extended periods of time, b) Trajectories were begun from 
identical configurations as in a, this time with noise added. 
With noise included, the metastability of the zero state is re- 
moved, and the gray region is now bounded away from zero. 



for D near the zero state. After a long time the sys- 
tem becomes polarized so that 5* <C 0, this allows us to 
adiabatically eliminate D+ to obtain 



iA S* + +m r S* z S + 
r-r i .-a mo* i 1Jz "r^A-^ ) 







(r + iA_)|5, 

■ 0(D 2 ) 



(11) 
(12) 



This linear stability analysis gives no conclusion about 
the stability of the zeros states. This result implies that 
within this model the stability of the zero state is only 
determined at higher order. This is a little surprising 



because at first glance Eq. 10 appears to have an attrac- 
tive force towards D z = 0. This arises from the same 
mechanism described in Ref. [22j however, a more careful 
treatment reveals that this effect actually cancels. Our 
simulations indicate that the nonlinear corrections make 
the zero state repulsive in the experimentally relevant 
parameter regimes. When we include the nuclear spin 
noise we shall show analytically that the system is re- 
pelled from the zero states. 
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FIG. 7: a) Asymptotic value of \D Z /S Z \ as a function of dot 
asymmetry with parameters chosen as in the location marked 
with an x in Fig.^] strongly in the instability regime. For each 
value of dot asymmetry R, we initialized fifty runs in a single 
initial spin configuration chosen from the thermal distribution 
(with D z = -0.72 and S z = -1.57). We plot the asymptotic 
value of D z /S z . The runs that ended with D z /S z greater 
(less) than shown are shown as red (blue) points. The circles 
(crosses) indicate the mean value of the red (blue) points, 
with error bars showing the standard deviations. The solid 
and dashed lines are given by Eq.|24|and Eq.[l3] respectively, 
b) As in (a), with parameters chosen in the location marked 
with an o in Fig. [4] strongly in the saturation regime. 



B. Effect of Nuclear Spin Noise 

1. Unequal Dots 

Our results that zero states are unstable to the growth 
of large difference fields in the presence of asymmetry in 
the size of the dots and nuclear noise can be be under- 
stood in the following heuristic picture first given in Ref. 
l23l We assume the nuclear spins have equal spin flip 
rates on the two dots, which is borne out by the analyt- 
ical and numerical calculations presented below. Then 
the build-up of the total Overhauser field S z is propor- 
tional to — (ge+g r ), where are the effective hyperfine 
interactions on the left (right) dot and the negative sign 
arise because nuclear spins are flipped down in the ex- 
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perimental cycles. Similarly D z grows as — {gt — g r ) so 
that the ratio 

D z /S z -+ (g t - g r )/{gi + g r ). (13) 

In this section we demonstrate a similar result within 
our full model. We assume homogeneous coupling and 
work in the high field, large J, limit where we can set 
Ao = A_ = in Prf. The local noise processes included 
in Eq. [7] give rise to a mean decay of the collective nu- 
clear spin variables and associated fluctuations , for 

L+(R+), defined by (T d (t)J^ t {tf)) n = 2Sl 2 d 5 dd ,5{t - t'). 
The semiclassical EOM for the nuclear spins reduce to 

L+ = g e r L Z (L + - R + )/2 - V L+ + y/tojPt, (14) 



Fo {L\ — Ri • Lj 



(15) 



and similarly for R, where r] is defined in Eq. [8j From 
Eq. 14 , we see that if we start in a zero state, Td will 
produce a fluctuation in D± , and the contribution to L z 
of the form —g^V^L\ results, in the long time limit, in 
L z <C — 1 and similarly for R z . Thus, \L Z /L Z \ <C 1 and 
we can treat L z , R z as static to find , {R\) n and 

(L_l • H±) n , which allow us to find the slow evolution of 
L z -> Rz- 

In particular, assuming L z , R z are constant we can 
write the closed set of equations for L + and R+ 



R + 



To 
2 



9iL z 
-g r R z 



~9fL z 
g r R z 



R + 



(16) 



Introducing the variables 



s + 



9eL z 
g r R z 
-1 



R + 



we find 

s+(t) = 

D+(t) = 



^1 f ^e-fo+^-'V/-^) 



(17) 

(18) 
(19) 



here 75 = —To{giL z + g r R z )/2 > 0. We can use this 
solution to calculate (£j_) n , (R\) n , and (L±-R±) n . 
For example to lowest order in 1/L Z , 1/R Z 



^/g 
' (1+p) 2 
gt + g r p 2 

2r] 



{91 + 9r)p 2 , 2p(g £ - g r p) 



275 



7s 



(20) 



where we have defined p = g^L z / g r R z1 g = [gg + g r )/2 
and used the fact that £l 2 d = g^jg in our units. 

Inserting this solution into the EOM for D Z1 S z gives 
reduced EOM for the slow, noise-averaged evolution of 



D z and S z . After some straightforward manipulations 
we arrive at 



S z 



gt v gt g r 
1g \S Z \ 2 



E 



S z 
D z 



where 5| = (gtL z +g r R z )/2 = -j s /^o and 

4R\(1- R)(l + R? -(1 + - R 2 ) 
and i? = grj gt- After rescaling time to 



Jo 



dt' 



9i V ge 9r 
9 \SI(t')\ 2 



(21) 



(22) 



(23) 



this becomes a purely linear system characterized by the 
matrix E. For all R > 0, this matrix has one positive 
and one negative eigenvalue; thus, it has one growing 
mode and one decaying mode. In the long time limit, 
both S z and D z will be proportional to their overlap with 
the growing mode. Thus D z /S z approaches a constant, 
which is easily found from E as 



1-R 2 



S z 2R+^J4R 2 + (1-R) 4 ' 



(24) 



In Fig. we compare this result and Eq. [l3]to the full 
numerics including all the parameters. The dot asym- 
metry favors D z /S z > 0, and the coherent instability 
mechanism does not prefer either sign. The competi- 
tion between the coherent and asymmetric mechanisms 
can be seen in the move toward zero of the points with 
D z /S z negative in Fig. [7^1; however, for larger asymme- 
tries (R > 0.5) all trajectories are seen to follow the di- 
rection of the natural asymmetry. Fig.[7]3 shows the same 
simulations performed in the saturation regime. As there 
is no coherent instability mechanism competing with the 
dot asymmetry, the sign of D z is determined by the asym- 
metry in all but the most symmetric dots. D z /S z is 
in good agreement with the simple prediction given by 



Eq. 13 and Eq. 24 



2. Identical Dots 

For identical dots the arguments given in the previous 
subsection break down; however, we shall now show that 
for certain parameters there still exists a mechanism for 
self-consistent growth of \D Z \. Growth of \D Z \ requires 
nonzero D±. For intermediate field and exchange, the 
A ,- contributions to become comparable to the r 
term. In particular, the AoD z z term acts as a source 
term for D± (see Eq. 25). Consequently, for weak enough 
noise D± will only be appreciable when \AoD z /TqS z \ is 
appreciable, which provides a self-consistency condition 
for the continued growth of D z . 

These properties of identical dots can be seen analyt- 
ically in the following limiting case: we assume a wave 



10 



function where the coupling takes two values, g\ ^> g 2 ,r\ 
and that initially —g2S z ^> gi \D Z \ ^> S± ~ 1 and 
D± ~ D z /S z <C 1. We denote the total angular momen- 
tum of nuclear spins in dot d with coupling constant g\~ 
by Jkd and assume J^ d ~ ~ J| d <C Jf d so that the 
majority of the polarization resides in the strongly cou- 
pled spins. We can write a closed set of equations for the 
evolution of D and S 

D + = gi ih_S z D + - gi iA D z S+ 

+g 2 5iA D z (J+ + J+)/2 - g 2 SiA_ D + (J^ + J| r )/2, 
5+ = -0iz(A o - A_)D Z D + + g 2 5i A D Z (J+ - J+)/2 

-g 2 5iA_D + (J* £ -J* r )/2 : 
D z = 9l Tm[A.D+S-] - g 2 Slm[A_D + (J- + J 2 ")/2], 
& = -^r ^l -^ 2 ^Im[A_D + (J 2 - - J 2 ")/2], 

= ±g2iA D z J+ d =f g 2 iA_D + J% d - rjJ+ d + / d , 
J$ d = ±g 2 Im[A_D + J-], 

where the top sign is for d = A_ = A_ — zTq, 5 = 
#i — #2, /d is a gaussian, white noise process derived 
analogously to J 1 ^ such that (fdfd) n ~ ^^cr 2 ! , and we have 
neglected to write the noise terms in the EOM for D+ 
and S+ because we have assumed they are higher order. 
Furthermore, we can neglect all terms proportional to 
g 2 D+J2 d because these are second order. This leads to 
the somewhat simpler set of equations 

D+ = gi iAS z D+ - gi iA D z S + (25) 

+ g 2 5iA D z (J+ + J+)/2, 

S + = - gi i{A -A_)D z D + (26) 

+ g 2 6iA D z (J+ - J 2 +)/2, 

j+ = ±g 2 iA D z J+ - riJ+ + / d , (27) 

Dz = #iIm[A_L> + #_] (28) 

S z = -giToD 2 ± , (29) 

These equations can be solved perturbatively in 
1/S Z ,1/D Z by the same method as in the previous sec- 
tion. The only difference in the structure of the two prob- 
lems is that in this case the source terms for D+ and S+ 
are proportional to J^ d instead of white noise; as a result 
we have to take into account the coherent evolution of 
the source term. We can expand the resulting EOM for 
D z in g\D z / g 2 S z to find the noise- averaged equation 

(rg + A2_-A A-W ^\ 3 (30) 

r§ + A2_ g 2 \\S z \) 

from which we see that the sign of Tq + A 2 _ — AqA_ 
determines whether or not there is continued growth of 
D z . Note that the perturbation theory breaks down as 
g 2 — >• 0. This reflects the importance of including the 




lo Sio m o 

FIG. 8: Phase diagram as in Fig.[4]except with noise added. 
The phase diagram is nearly identical. See Table I for param- 
eters. 

coherent evolution of in solving for the dynamics. 
Without g 2l we would have found D z = 0. This phase 
boundary is shown as the dashed line in Fig. [3] In Fig. [8] 
we show the phase diagram as a function of cycle fre- 
quency and inverse magnetic field, where we see qualita- 
tively the same behavior as Fig. [4] 



V. RELEVANCE TO OTHER CENTRAL SPIN 
SYSTEMS 

Although this work has focused on lateral double quan- 
tum dots in GaAs, the methods, and some of the results, 
can be ap plied to vertical double dotd 3 ^, InAs quantum 
dot d 35 * 36 *, silicon based quantum dots^, and NV-centers 
in diamond^. A few important differences for these other 
central spin systems are that the sign of the electron g- 
factor may be positive (compared to its negative sign in 
GaAs) and the spin-orbit coupling can be much larger 
in other systems than it is in GaAs 43 . The results pre- 
sented in the paper are not dependent on the sign of the 
^-factor. Changing the sign would reverse the direction 
of the nuclear polarization from negative to positive, but 
all of our analysis would carry through essentially un- 
changed. The competition between spin-orbit coupling 
and DNP is more dramatic and can have a qualitative 
effect on the polarization dynamics for large spin-orbit 
coupling^. 



VI. CONCLUSIONS 

We have shown that dynamic nuclear polarization ex- 
periments in double quantum dots give rise to a rich set 
of phenomena. We have developed novel numerical and 
analytical methods to theoretically describe such dynam- 
ics. We find that the competition between polarization, 
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noise processes and coherent evolution mediated by the 
electrons allows one to carefully control the final nuclear 
spin state in the two dots. The main results of the paper 
that the dynamics are dominated by either rapid sat- 
uration of polarization or an instability to the growth 
of large difference fields are consistent with the obser- 
vations of the experiments reported in Refs. [18j [13] and 
124} however, we see evidence that the dynamics is much 
richer as the experiments have not resolved whether or 
not the instability to large difference fields results from 
dot asymmetry or coherent electron- nuclear interactions. 
These two cases could be experimentally distinguished 
by measuring the sign of D z in a given double dot. In 
addition, we see that the zero states may be experimen- 
tally observable as metastable states in certain parameter 
regimes, indicating that there is still much to explore in 
the polarization dynamics of double quantum dots. 
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Appendix A: Parameters Used in Simulations 

In table I below we provide a summary of the param- 
eters used in the simulations for each figure. 



Appendix B: $ variables 

In this appendix we describe a systematic approach to 
coarse graining the electron wavefunction in solving the 
semiclassical equations of motion, which we refer to as the 
Independent Random Variable Annular Approximation 
(IRVAA). We construct a sequence of discretizations of 
the wavefunction for which we can provide a rigorous 
bound on the error in time evolution compared to the 
exact solution. In the process we also introduce a new set 
of statistically independent nuclear spin variables, which 
are a convenient basis for numerical simulations. 

We see from Eqs. [4] and [6] that the semiclassical evolu- 
tion of each spin depends only on the vectors L and R 
(or equivalently on D and S). That is, if we know Pd(t) 
(which depends only on L and R), then we can solve for 
the dynamics of the entire system. However, even if we 
know Pd(t), if we look at the equation of motion for L we 
find that it generates an infinite hierarchy of equations 



dL 

~dt 



where we defined L* = ^ k g^iki- Now L* couples to the 
variable J2 k glihi and so on. 

To find an approximate solution to the dynamics we 
would like to find an effective method to truncate this 
infinite hierarchy of equations. For simplicity we focus 
on the case where Pi is only a function of L, reducing 
it to a single dot problem, and drop the dot indices in 
the following discussion. We also work in the continuum 
limit, which is defined by a nuclear angular momentum 
density I(r,t) = J2k*k(t)5(r - r k ). 

Each v ariab le in the hierarchy of equations of motion 
(as in Eq. |B1|) can be expressed as an integral 



<S>(t) = j d d rg(r)^g(r))I(r. t), 



(B2) 



where (p(x) is a polynomial in x. That is, there is a one- 
to-one correspondence between polynomials cp(x) and the 
variables in the EOM. For example, L corresponds to 
<t>{x) = 1. 

We would like to think of a truncation procedure as 
any procedure that provides a reduced, self-consistent set 
of equations describing the evolution of P, equivalently 
L. We make a formal definition of a truncation pro- 
cedure as a procedure producing a set of variables 
k = 1, . . . , M, of the form above and an M x M matrix 
Q, such that 3>i = L and 



dt 



Since we always constrain $1 = L, we always have 
<Mz) = l. 

To construct a convenient basis of nuclear spin vari- 
ables we first define a norm (•) based on the statistical 
average of a nuclear spin variable in the infinite temper- 
ature ensemble, i.e. 



d d rd d r'g 2 (r) V {g{r))^{g(r')) (I(r) ■ I(r% 



1(1+1) 



J d d rg 2 (rMg(r))iP(g(r)) 



(B3) 



= Pi x L* 



(Bf) 



where a is the lattice spacing, (-) e is the ensemble average 
over the initial thermal state and we took (J(r) • I(r')) = 
1(1 + l)S(r — r')/a d . Now we can construct an or- 
thogonal set of polynomials with respect to this norm 
by using the standard Gram- Schmidt procedure start- 
ing from the polynomial 1. This gives a set of orthog- 
onal polynomials (p k and associated nuclear spin vari- 
ables <3>/c = J d d r g(r)(pk(g(r))I(r,t), which are statis- 
tically independent in the infinite temperature ensemble 
(i.e., (<$> k ' *&i) = 3f^<5fcz)and satisfy <I>i = L. 

The equations of motion (EOM) for these variables can 
be written as 

$ n = PX Q nm ^m (B4) 

where the matrix Q mn is a tridiagonal matrix defined by 
the recurrence relations 

X(p n (x) = Qnn-Wn-1 + Qnn^n + Qnn+Wn+l (B5) 
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TABLE I: Parameters used in the simulations shown in the figures of this paper. 
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and we used the fact that xip n (x) only has a non-zero 
overlap with <p n and ip n ±i. 

We now define an M th order truncation procedure with 
respect to the variables by setting Qmm+i = 0. The 
central result of this appendix is encapsulated by the fol- 
lowing theorem for this truncation procedure. 

Theorem: For a given wavef unction g(r) and e > 0, 
the above truncation procedure at order M will produce 
an effective L M (t)such that \L(t) - L M (t)\ < e for all 
t < tM, where tu is a time scale that increases linearly 
with M and L(t) is the exact result for the untruncated 
system. 

We begin our analysis by proving that any truncation 
procedure is equivalent to a discretization of the function 
g(r) (i.e., an annular approximation), by which we mean 
a representation of L as 



M 



L = ^2g(r k )I k 



(B6) 



k=i 



where I k is a rescaled nuclear spin variable associated 
with position r&. 

The reverse implication is clear because if we start with 
such a discrete representation, then the variable associ- 
ated with the polynomial 



A I 



w ( x ) = ni x -g( r k)} 



k=l 



is identically zero. That is, if there are only M discrete 
spins in the system, then there are only M statistically 
independent variables <f>k in the system, and <E>m+i is 



naturally zero. This result naturally truncates Eq. B4 



degree less than M and its associated set of spin vari- 
ables, then we can obtain a finite, self-consistent set of 
equations for the evolution of L. 

The forward implication follows along similar lines. If 
M — 1 is the maximal degree of the set of polynomials 
{cpk(x)} associated with the truncation variables {&k} 
and $m is the spin variable corresponding to this poly- 
nomial, then, when we compare to the continuum limit, 
we find that the statement that d&M/dt does not couple 
to higher degree polynomial variables implies the exis- 
tence of a degree-M polynomial w(x) such that 



/ 



d d rg(r)w(g(r))I(r,t) = 0, 



for any J(r, t). The existence of such a polynomial imme- 
diately implies that we can represent L in the discretized 



form of Eq. B6 



We have now reduced the problem of finding an opti- 
mal truncation procedure to the problem of finding an 
optimal discretization procedure for integrals of the form 



/ 



d d rg(rMg(r))I(r,t), 



Consequently, if we consider any basis of polynomials of 



where <p(x) is a polynomial in x. Fortunately, this 
last problem is solved through the theory of Gaussian 
quadrature.^ First, though, we assume that our func- 
tion g(r) is spherically symmetric so that we can write 
our integrals as effective one-dimensional integrals with 
respect to the rescaled angular momentum density 

J(r, t) = J dtt a d ~ x N(r) I(r, ft, t) /S(d) (B7) 

where ft parameterizes the surface of a d-dimensional 
sphere, a is the lattice spacing, S(d) is the surface 
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area of a unit sphere in d dimensions, and N(r) = 
S{d)r d ~ 1 / a d ~ x is the number of nuclear spins at radius 
r; for example, in two dimensions N(r) = 2itr/a. The 
ensemble average of I(r,t) is given by (I(r) - I(r')} = 
I{I+l)N{r)5{r-r')/a. 

To begin constructing our Gaussian quadrature rules 
we rewrite 



*(*) = / drN(r)g 2 (rMg(r)) 
Jo 



I(r,t) 

N(r)g(r) 

o N(g 1 (x))x 



(B8) 



where x = g(r) and u{x) = ^ \g~ 1 (x) (x))x 2 is the 
weight function. Standard results in the theory of numer- 
ical integration imply the existence of a set of orthogonal 
polynomials, ip n , with respect to the inner product 



(/ 



h)= f 

Jo 



dxoo(x) f(x) h(x) 



(B9) 



such that, for any function f(x), the M th order quadra- 
ture approximation is given by 



Jo 



M 



dxu(x) f(x) « y^Ukfjxk), (BIO) 



k=i 



where x k are the zeros of cfM and the weights uo k are de- 



termined by the condition that Eq. BIO is exact for all 
polynomials of degree strictly less than 2M. The error in 
this formula decreases exponentially in M, or better, pro- 
vided that / is smooth™ In addition, these polynomials 
are exactly the ones we used to construct our truncation 
procedure. Consequently, our truncation procedure de- 
fined above is equ ivalent to approximating L in quadra- 
ture as in Eq. 



B6 



with I k = w k I(r k ,t)/glN(r k ). 
To prove the theorem we first note that from the 
definition \P(L)\ < 1 for all L. Now let p > be 
such that \P(L) - P{L')\< p\L - L'\ for all L and 
L . We define L n (t) = J d d rg n (r)I(r,t) and L™(t) 
is the solution for the equivalent variable in the trun- 
cated system of equations. To provide bounds on the 
error propagation we define S^f (t) = \L n (t) — L^f (t)\. 
We work in time units where max r g(r) = 1 and let 
b = max n5 t|L n (t)|< J d d rg(r)(I + 1). Now it is straight- 
forward to show that 



5?<pb5f 



M 



(l+P<5f)C+i<CW 



M 



C+i) 



(Bll) 



where £ = max(p6, 1 + pe) and, by assumption, we are 
restricted to short enough times that S^ 1 < e. By con- 
struction, 6}f(0) = for n < M while for n > M 
dff is bounded by the quadrature error on the integral 
J d d rg n (r)I(r,0), which is less than ce~ M for a con- 
stant c independent of M. Using Eq. Bll| we can then 
bound the error on Sf 1 < ce~ M (e 2( ^ t — 1). This im- 
plies that the time to make an error of size e scales as 
(l/2C)log(£:e M /c + 1) ~ (M - \ogc/e)/2( for large M. 
This proves the theorem. 



For the two dimensional Gaussian g(r) oc e _r2//2cr2 the 
weight function w(x) = x and the associated orthogonal 
polynomials are the Jacobi polynomials. The matrix Q 
is then given by standard recurrence relations for Jacobi 
polynomials. Once the recurrence relations are known, 
one can work with the ^-variables without converting 
between the original nuclear spin variables because the 

variables were defined such that they are initially sta- 
tistically independent. This is a convenient numerical 
approach for these types of central spin problems, and it 
was used in all of the numerics in this work. 



Appendix C: Multiple Nuclear Species 

In this appendix we include the effects of multiple nu- 
clear species in our simulations and find that the main 
results for both asymmetric and identical results carry 
through much the same. First we show how to include 
multiple species in terms of the collective ^-variables and 
then we present the simulation results. 

When multiple species are taken into account we must 
include the Larmor precession of the nuclear spins. In 
this case the EOM take the form 



L kd 



: j e b a v \i>kd\ Pd x /; 



kd 



kdi 



(CI) 



where a is a species index, u a = j a B ext T/r a is the ef- 
fective Larmor frequency, b a is the bare hyperfine field 
of species a, j a is the gyromagnetic ratio of species a, 
B ext is the external magnetic field, and we have explicitly 
included the factor T/r a , where T is the total time of the 
nuclear pump cycle and r a is the adiabatic sweep time. 



TABLE II: Relative population of the nuclear species x a , ef- 
fective hyperfine field due to species ol 0a. , and the gyromag- 
netic ratio 7c, for the three nuclear species in GaAs. 





75 As 


69 Ga 


71 Ga 




1 


0.6 


0.4 


ba (T) 


-1.84 


-1.52 


-1.95 


sy ( kH^A 

la y mT j 


45.96 


64.39 


81.81 



We introduce the projector function 7r kd , such that 
ir kd = 1 if there is species a in unit cell k and oth- 
erwise. This allows us to write 



L = ^2leb a V0 \Ai\ 2 ^k£ J k£ 
k,a 



g M t^m 1 ki- 



rn 



VJ2a ba X a ka 

Here we have defined Qp to be the standard deviation of 
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t 



> D z /S z positive 
< D z /S z negative 



g r /gi 




A_/A 



FIG. 9: a) As in Fig. [7] with parameters chosen as in Fig. 2 of 
Ref. 1231 except with three species. Due to the computational 
cost of running three species of spins, simulations were run 
for only 10% as long, and the range of D z /S z is larger as a 
result. The trend that D z /S z is in good agreement with the 
single-species prediction is clearly visible, b) Phase diagram 
with multiple species and mo = 0. 



L M in the infinite temperature state, explicitly 

<J& ■ I%', d ,) = 1(1 + 1)<W<W<W, (C3) 
^<L 2 )/3 = ]T 7e 2 ^ Q t, 2 |^ (C4) 



where x a = (7r% d ) is the relative proportion of species a 
on the sites it can occupy, gkd vq \ipkd\ 2 are chosen to 
satisfy ^2 k g\ 1(1 + 1) = 3, and / is the total spin of a 
single nuclear spin (7 = 3/2 for all species in GaAs). 
We define the variables 



n = ^=Y,9ki^n{9ki)^I^ (C5) 



where (Pni x ) are defined as in Appendix B and are inde- 
pendent of the species, i.e. ¥>q(x) = 1 and 

1>L V&(9kd) vLiShd) HI + 1) = 3 S nm . (C6) 



These definitions have the implication that 
= <W<W<W, and we can draw ini- 
tial values for each of them from a normal distribution. 
Furthermore, we can express 



ft/ 



(C7) 



All these definitions are equivalent for the right dot. 
In these variables the EOM take the form 



+ e n+1 $« +1 ) -to a z x *° 



(C8) 



where we have used the definition TV -1 = max^ \^kd\ 2 
to represent the number of nuclear spins with which 
the electron has significant overlap. For a two di- 
mensional gaussian wave function we have N = 

2/3£ a ^7eW+l)/" 2 

In Fig. [9] we include the three nuclear species in the 
simulation and show that qualitatively the results from 
the single species case still hold. Fig. [9^ shows the asymp- 
totic ratio of D z /S z as the relative dot sizes are varied, 
where we see good agreement with the simple prediction 
given in the introduction. In Fig.[9]3 we extract the phase 
diagram in the simplified model with only Ao,_ and To 
non-zero, as in the model of Ref. 23. As in the single- 
spin case, we find a saturation regime at high values of 
r /A and an instability regime at lower values. Unlike 
in the single-spin case, the saturation regime does not 
broaden at higher values of A_/Aq. The dashed line is 
the same as that in Fig. [3J showing the simple predic- 
tion for the phase boundary with a single species, from 
Ref. |23j The lower-left side of the phase diagram (the re- 
gion most easily reached in experiments) is well-described 
by this prediction, even with multiple species. 
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